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Abstract 

This  paper  addresses  the  development  of  a  temperature-dependent  constitutive  model  for  relaxor  ferroelectrics 
which  is  based  on  the  assumption  that  the  material  is  comprised  of  an  aggregate  of  micropolar  regions  having  a 
range  of  Curie  temperatures.  The  diffuse  transition  behavior  of  the  material  is  due  to  its  chemical  heterogeneity, 
and  thermodynamic  models  for  the  microregions  are  developed  by  considering  near  neighbor  interactions  for  varying 
cation  ratios.  The  resulting  micropolar  model  can  be  used  to  predict  the  saturation  polarization  and  distribution  of 
regions  as  a  function  of  temperature.  Hysteresis  below  the  freezing  point  is  incorporated  through  the  quantification 
of  energy  required  to  bend  and  translate  domain  walls  pinned  at  inclusions  inherent  to  the  material.  The  resulting 
ODE  model  quantifies  the  constitutive  nonlinearities  and  hysteresis  exhibited  by  the  materials  through  a  wide  range 
of  temperatures  and  input  drive  levels.  The  predictive  capabilities  of  the  model  are  illustrated  through  a  comparison 
with  PMN-PT-BT  data  collected  at  temperatures  ranging  from  263°  K  to  313°  K. 
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1.  Introduction 

Actuators  constructed  from  the  relaxor  ferroelectric  material  Pb(Mg1//3,  Nb2/3)03  (PMN)  provide  higher  induced 
displacement  and  lower  hysteresis  than  analogous  devices  constructed  from  conventional  piezoelectrics.  When  used 
as  drivers  for  underwater  sonar  projectors  [12,  15],  PMN  has  demonstrated  at  least  a  6  dB  increase  in  acoustic  source 
level  over  Pb(Zr,  TijCR  (PZT).  When  employed  in  multi-layered  arrays  for  deformable  mirrors  [6,  9,  10],  PMN  has 
increased  mirror  sensitivity  from  0.3  nm/V  to  12.0  nm/V  relative  to  monolithic  PZT. 

Whereas  PMN  has  enabled  improved  system  performance,  its  nonlinear  electro-mechanical  behavior  has  also 
added  more  complexity  to  the  design  process.  PMN  is  not  poled  during  the  manufacturing  process  because  its  induced 
strain  is  electrostrictive.  A  designer  must  select  an  electric  field  bias  so  the  device  operates  on  one  wing  of  the  strain- 
electric  field  curve.  Additionally,  PMN’s  induced  strain  and  hysteresis  are  highly  dependent  on  both  temperature  and 
applied  frequency.  When  driven  at  high  fields  and  high  frequency,  PMN  will  heat  due  to  dielectric  hysteresis.  Since 
PMN’s  hysteresis  decreases  with  temperature,  actuators  reach  an  equilibrium  temperature  when  heat  dissipation  due 
to  conduction,  convection  and  radiation  balances  the  internal  heat  generated  by  the  ceramic.  The  final  equilibrium 
temperature  depends  on  the  configuration  of  the  actuator,  its  electric  drive  cycle,  and  its  surrounding  environment. 
Shankar  and  Horn  [16]  predicted  a  40°  C  rise  for  a  flextensional  sonar  transducer  submersed  in  water  using  a  simple 
heat  generation  model  for  PMN.  These  predictions  agreed  quantitatively  with  experimental  measurements  made  on 
an  actual  transducer. 

In  order  to  design  PMN  compositions  that  minimize  this  effect,  we  need  temperature-dependent  constitutive  laws 
to  use  as  transducer  design  tools  for  simulating  the  self-heating  phenomenon.  In  recent  years,  nonlinear  constitutive 
models  have  been  developed  for  PMN  to  simulate  their  response  and  internal  stress  states  at  a  fixed  temperature. 
Horn  and  Shankar  [11]  proposed  an  isothermal,  electrostriction  model  that  includes  stress  effects  on  the  dielectric 
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response  and  nonlinear  saturation  of  the  polarization  at  high  fields.  Subsequent  experiments  by  Brown  et  al.  [2] 
demonstrated  that  this  model  accurately  predicts  the  coupled  electro-mechanical  behavior  of  PMN  under  isothermal 
conditions  in  the  absence  of  hysteresis.  Smith  and  Horn  [13]  added  hysteresis  to  this  constitutive  law  by  modeling 
domain  motion  and  pinning.  This  new  model  simulated  the  electrical  behavior  at  temperatures  significantly  below 
the  Curie  temperature  when  large  domain  structures  have  formed.  In  a  parallel  effort,  Glazounov  et  al.  [7]  used 
the  same  polarization  model  as  Horn  and  Shankar  to  analyze  the  temperature  response  of  PMN.  They  assumed 
that  the  material  consists  of  micropolar  regions  with  a  distribution  of  local  transition  temperatures.  They  deduced 
the  saturation  polarization  of  the  aggregate,  the  density  of  the  micro  regions  and  the  distribution  of  the  transition 
temperatures  from  experimental  data  over  a  broad  temperature  range.  However,  their  model  did  not  account  for 
interaction  between  the  micropolar  regions. 

In  this  paper,  we  extend  the  work  of  Glazounov  et  al.  to  develop  a  static,  temperature-dependent  constitutive 
model  for  relaxor  ferroelectric  materials.  Like  Glazounov  et  al.,  we  assume  that  the  relaxor  material  consists  of 
superparaelectric,  micropolar  regions  with  a  diffuse  spectrum  of  Curie  temperatures.  However,  we  use  a  simple  Ising 
model  with  near  neighbor  ion  interaction  to  represent  the  thermodynamics  of  the  individual  micropolar  regions. 
Based  on  the  “random-layer”  model  of  spatial  B-site  ordering  for  relaxor  ferroelectrics,  we  calculate  the  dispersion 
of  local  Curie  temperatures.  The  micropolar  model  predicts  the  saturation  polarization  and  the  density  of  the  polar 
regions  as  a  function  of  temperature.  A  macroscopic  model,  which  includes  interaction  between  the  micropolar 
regions,  was  constructed  with  these  parameters  to  simulate  polarization  response  of  the  relaxor.  Finally,  a  domain 
wall  motion  model  based  on  that  proposed  by  Smith  and  Horn  [9]  was  added  to  incorporate  hysteresis  at  temperatures 
significantly  below  the  macroscopic  Curie  temperature.  The  model  predictions  agree  qualitatively  with  experimental 
measurements  on  PMN  through  a  large  temperature  range. 

1.1.  Physical  Mechanisms  for  Relaxor  Behavior 

“Normal”  ferroelectrics  exhibit  a  sharp  ferroelectric-to-paraelectric  phase  transition  and  piezoelectrics  like  PZT 
operate  at  temperature  well  below  this  transition.  In  contrast,  relaxor  materials  possess  a  diffuse  phase  transition 
over  a  broad  temperature  range  [7,  20].  Viehland  et  al.  [21,  22]  have  demonstrated  that  relaxors  deviate  from  the 
Curie- Weiss  law  at  temperatures  above  the  temperature  of  dielectric  maximum,  Tm,  and  that  they  exhibit  non- 
Debye  relaxation  cumulating  in  polarization  freezing  at  a  temperature  below  Tm .  Since  relaxors  typically  operate  at 
temperatures  within  the  transition,  the  material  behavior  strongly  depends  on  the  temperature  and  frequency  of  the 
applied  electric  load  as  illustrated  in  Figure  1.  The  transition  temperature  is  often  tailored  to  a  specific  application 
by  mixing  PMN  with  lead  titanate  (PT)  and  additional  minor  dopants,  such  as  barium  titanate  (BT)  and  strontium 
titanate  (ST).  The  PMN  family  has  the  perovskite  structure  ABO3  shown  in  Figure  2.  The  A  sites  are  lead  anions, 
while  the  B  sites  are  a  mixture  of  low  valence  cations  and  high  valence  cations  (1/3  Mg+2  and  2/3  Nb+5  in  the  case 
of  PMN).  The  spatial  ordering  of  the  multiple  B-site  cations  leads  to  the  diffuse  transition  of  a  relaxor.  Smolensky 


Figure  1.  Dielectric  response  of  a  PMN-PT-BT  (10%  PT,  3%  BT)  at  frequencies  ranging  from  0.1  kHz  to  100  Hz. 
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Figure  2.  Atomic  structure  of  a  perovskite  relaxor  ferroelectric  A(B’,B”)03.  The  B’  and  B”  sites  indicate  the  cation 
ordering  for  the  “random  layer”  model. 


postulated  a  dispersion  of  cations  which  creates  Mg:Nb  rich  ferroelectric  regions  surrounded  by  pure  Nb  sites  [20]. 
Each  region  has  a  local  phase  transition,  so  the  aggregate  possesses  a  spectrum  of  phase  transitions.  Cross  [5] 
later  proposed  the  more  complete  picture  of  relaxor  behavior  shown  in  Figure  3.  He  postulated  that  nanoscale  size 
(3-5  nm)  polar  regions  form  with  decreasing  temperature.  These  small  regions  have  a  low  energy  barrier  so  they 
fluctuate  with  thermal  agitation  in  a  manner  analogous  to  superparamagnetism.  Each  region’s  dipole  moment  has 
random  orientation;  the  macroscopic  polarization  is  zero  but  the  system’s  RMS  dipole  moment  is  nonzero.  Initially, 
the  relaxor  obeys  the  Curie- Weiss  relationship  with  a  very  high  Curie  temperature  (398°  K  for  PMN).  With  further 
cooling,  the  thermal  fluctuations  slow  down  and  interaction  between  regions  freezes  the  polarization  in  a  manner 
similar  to  magnetic  spin  glasses.  The  material  deviates  from  the  Curie- Weiss  relationship  and  possesses  a  strong 
frequency  dependence.  Eventually,  the  material  develops  a  permanent  remanent  polarization  and  large  domains 
develop  as  in  a  “normal”  ferroelectric. 

Clearly,  the  ordering  of  the  B-site  cations  defines  the  diffuse  nature  of  the  transition  and  two  cation  ordering 
models  have  been  proposed.  Diffraction  studies  indicate  that  the  polar  regions  have  a  1-to-l  distribution  of  two 
different  cation  sites  arranged  in  a  face  center  lattice.  The  two  sites  are  denoted  as  B’  and  B”  as  depicted  in 
Figure  2.  In  the  “space  charge”  model  for  PMN  [8,  22],  Mg  exclusively  occupies  the  B’  sites  in  the  ordered  regions, 
while  Nb  exclusively  occupies  the  B”  sites.  The  ordered  regions  thus  carry  a  net  negative  charge  which  should  inhibit 
domain  growth.  In  the  “random  layer”  model  [1,  3],  Mg  and  Nb  in  a  2-to-l  ratio  randomly  populate  the  B’  sites, 
while  Nb  exclusively  occupies  the  B”  sites.  In  this  model,  the  ordered  regions  are  charge-balanced  so  domains  should 
coarsen  during  heat  treatment.  Recent  experiments  by  Akbas  and  Davies  [1]  have  shown  that  domain  growth  does 
occur  during  annealing  in  Pb(Mg1(/3  Ta2/3)03  (PMT),  which  supports  the  “random-layer”  model  of  the  B-site  cation 
distribution. 


2.  Isothermal  Hysteresis  Model 

The  methodology  for  constructing  a  temperature-dependent  model  for  hysteresis  and  constitutive  nonlinearities 
in  relaxor  ferroelectrics  is  motivated  by  that  employed  in  [17,  18]  for  fixed  temperatures.  The  first  step  comprises 
the  development  of  an  appropriate  anhysteretic  model  which  can  be  single  or  multi-valued  according  to  the  choice 
of  parameters.  While  the  multi-valued  anhysteretic  curve  incorporates  a  form  of  hysteresis,  the  transition  between 
positive  and  negative  remanence  values  is  typically  steeper  than  that  observed  in  polycrystalline  materials  at  tem¬ 
peratures  approaching  the  Curie  point.  The  second  component  of  the  model  quantifies  reversible  and  irreversible 
domain  wall  losses  due  to  bending  and  translation  of  domain  walls  pinned  at  inclusions  inherent  to  the  material. 
The  combined  model  characterizes  both  steep  and  gradual  transitions  between  remanence  values  and  hence  provides 
predictive  capabilities  through  a  range  of  operating  conditions. 

We  summarize  here  the  model  developed  in  [17,  18]  to  illustrate  the  methodology  and  to  indicate  the  extensions 
required  to  accommodate  changing  temperatures.  An  important  attribute  of  the  model  is  the  property  that  model 
parameters  can  be  asymptotically  correlated  with  physical  properties  of  the  data  including  the  coercive  field,  rema¬ 
nence  point,  and  the  differential  susceptibility  of  the  data  at  these  points.  Certain  asymptotic  relations  developed  in 
[19]  are  summarized  to  highlight  this  attribute  and  motivate  components  of  the  temperature-dependent  model. 
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Increasing  Temperature 


Ferroelectric:  Micropolar  re¬ 
gions  interact  and  material  ex¬ 
hibits  remanent  polarization. 


Superparaelectric:  Micropolar 
regions  develop  but  have  random 
orientation  so  no  remanence. 


Paraelectric:  Micro-regions  are 
nonpolar  and  material  exhibits 
zero  macroscopic  polarization. 


Figure  3.  Polar  mechanisms  in  relaxor  ferroelectric  materials  as  postulated  by  Cross  [5].  In  the  superparaelectric 
phase,  the  exhibits  a  nonzero  RMS  dipole  moment  but  has  zero  macroscopic  polarization. 


2.1.  Hysteresis  Model  for  Fixed  Temperatures 


The  fixed-temperature  hysteresis  model  is  constructed  in  two  steps;  (1)  the  characterization  of  the  anhysteretic 
(hysteresis-free)  relation  between  the  input  field  E  and  polarization  P,  and  (2)  the  quantification  of  energy  required 
to  bend  and  translate  domain  walls  pinned  at  inclusions  in  the  material. 

One  approach  used  to  model  the  anhysteretic  polarization  Pan  is  to  employ  Boltzmann  statistics  to  quantify  the 
probability  that  dipoles  occupy  certain  energy  states.  As  detailed  in  [17,  18],  the  balance  of  electrostatic  and  thermal 
energy  with  varying  assumptions  regarding  possible  moment  orientations  yields  the  Ising  spin  and  Langevin  relations 


Pan  =  Ps  tanh 


P  =  P 

L  an  —  x  s 


(i) 


TV’  rr l  17  rr i 

where  a  =  —  in  the  Ising  spin  relation  and  a  =  the  Langevin  model.  The  parameters  PS,TC  and  Eq  respectively 

denote  the  saturation  polarization,  the  Curie  temperature  for  the  material,  and  a  scaling  electric  field.  The  effective 
field  at  the  domain  level  is  taken  to  be 


Ee  —  -h  cxP  +  Pa  . 


(2) 


The  component  aP  incorporates  the  effects  of  interdomain  coupling  while  Pa  quantifies  the  electromechanical  cou¬ 
pling  due  to  applied  stresses.  For  fixed-temperature  material  characterization,  the  parameters  a,  a  and  Ps  are 
estimated  using  either  the  asymptotic  relations  developed  in  [19]  or  a  least  squares  fit  to  data. 

The  second  component  of  the  model  quantifies  the  energy  required  to  translate  and  bend  domain  walls  pinned 
at  inclusions  in  the  material.  As  detailed  in  [17,  18],  the  irreversible  changes  in  polarization  due  to  domain  wall 
translation  are  incorporated  through  the  quantification  of  the  energy  required  to  reorient  dipole  moments  thus 
yielding  the  relation 


dPi7 


dE 


=  8 


kS  Oi  ( Pan  Pirr ) 


(3) 


for  the  irreversible  polarization.  The  parameter  6  =  sign(Fi)  ensure  that  the  energy  required  to  break  pinning  sites 
always  opposes  changes  in  the  polarization.  The  physical  observation  that  polarization  changes  after  a  reversal  in 
field  direction  are  reversible  motivates  the  incorporation  of  the  parameter  6  which  is  defined  to  1  if  {dE  <  0  and  P  > 
Pan}  or  {dE  <  0  and  P  >  Pan}  and  is  0  otherwise.  The  parameter  k  is  defined  by  k  =  where  n  denotes  the 
average  density  of  pinning  sites,  p  is  an  average  dipole  moment,  and  (£*•)  is  the  average  energy  for  180°  domain  walls. 
While  k  was  considered  constant  in  the  fixed-temperature  analysis  of  [17,  18,  19],  it  is  clearly  temperature-dependent 
due  to  the  dependence  of  {£n)  on  the  thermal  energy. 

The  second  component  of  the  polarization  is  the  reversible  polarization  which  models  the  effects  of  domain  wall 
bending.  To  a  first  approximation,  this  is  modeled  by  the  relation 


Prev  PP'vn  Pit 


(4) 
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where  c  is  a  parameter  which  is  estimated  for  a  given  material  [17,  18,  19]. 

The  total  polarization  is  given  by 

P  —  Pirr  +  Prev  ■  (5) 

To  implement  the  model,  the  effective  field  for  a  given  field  and  irreversible  polarization  level  is  computed  using 
(2).  This  effective  field  value  is  then  employed  in  (1)  to  compute  the  corresponding  anhysteretic  polarization.  The 
subsequent  irreversible  polarization  is  determined  by  numerically  integrating  (3).  The  total  polarization  is  then 
specified  by  (5). 


2.2.  Asymptotic  Parameter  Relations 


The  previous  model  requires  the  parameters  a,a,c,Pg  and  k.  Algorithms  for  determining  initial  estimates  of 
these  parameters  using  measured  values  of  the  remanence  polarization,  coercive  field,  peak  polarization  and  field, 
and  the  differential  susceptibilities  at  these  points  were  developed  in  [19].  We  summarize  here  those  asymptotic 
relations  which  pertain  to  the  development  of  the  temperature-dependent  model. 

Let  Xan  and  Xc  respectively  denote  the  differential  susceptibilities  of  the  anhysteretic  polarization  at  E  =  0  and 
the  total  polarization  at  E  =  Ec.  The  slopes  of  the  polarization  curve  before  and  after  tip  reversal  are  denoted  by 
X^  and  Xm  ■  As  developed  in  [19],  estimates  of  the  parameters  a  and  k  are  given  by  the  asymptotic  relations 


(i) 


Ps 

^ Xan 


(ii)  k  m  Ec  (For  soft  materials)  . 


The  saturation  polarization  Pg  can  quite  accurately  be  estimated  directly  from  high  drive  level  data  while  the  slope 
Xan  is  approximated  by  the  slope  at  coercivity  Xan  ~  Xc-  Details  concerning  these  relations  as  well  as  the  complete 
algorithm  for  estimating  the  full  parameter  set  are  provided  in  [19]. 


3.  Temperature-Dependent  Hysteresis  Model 

From  a  physical  perspective,  the  parameters  a,a,Pg,k  and  c  are  temperature-dependent  so  the  model  (5)  is 
valid  only  in  operating  regimes  in  which  temperatures  exhibit  minimal  variation.  In  this  section  we  address  the 
temperature-dependence  through  the  consideration  of  the  relaxor  ferroelectric  as  an  aggregate  of  micropolar  regions 
having  differing  densities  of  Mg  cations  and  hence  differing  ferroelectric  properties.  The  development  of  the  anhys¬ 
teretic  relations  is  motivated  by  Glazounov  et  al.  [7]  and  follows  the  theory  developed  by  Horn  and  Shankar  [13]. 
The  temperature-dependence  in  the  pinning  losses  is  incorporated  through  energy  arguments  analogous  to  those 
employed  in  the  theory  of  dislocation  plasticity  [14].  The  combined  model  can  be  used  to  characterize  hysteresis  in 
relaxor  ferroelectrics  throughout  the  range  of  temperatures  encountered  in  operation. 

3.1.  Anhysteretic  Polarization 

As  noted  previously,  relaxor  ferroelectrics  are  assumed  to  consist  of  micropolar  regions  whose  transition  temper¬ 
atures  and  ferroelectric  properties  are  determined  by  variations  in  the  concentration  of  cations  (Mg2+  and  Nb5+ 
in  PMN)  which  comprise  the  material.  This  component  of  the  model  is  constructed  in  two  steps.  In  the  first,  we 
consider  the  thermodynamics  of  the  micropolar  regions  in  order  to  predict  the  saturation  polarization  and  density 
of  Mg  cations  as  a  function  of  temperature.  Secondly,  the  interaction  between  regimes  is  quantified  to  construct  a 
macroscopic  model  which  predicts  the  dielectric  properties  and  polarization  of  the  material. 

For  modeling  purposes,  we  make  the  following  assumptions  which  are  depicted  in  Figure  4. 

(i)  Each  region  contains  a  total  of  No  Mg  cations. 

(ii)  The  density  of  Mg  cations  X  is  uniform  throughout  each  region  but  varies  between  regions.  As  a  result  of 
Assumption  (i),  the  volume  V)  is  each  region  is  inversely  proportional  to  Xt. 

(iii)  The  strength  of  interactions  between  regions  is  directly  proportional  to  the  density  of  nearest  neighbors. 

(iv)  Below  its  local  transition  temperature,  each  region  acts  as  a  single  dipole  and  exhibits  a  remanent  polarization. 
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Figure  4.  Relaxor  ferroelectric  consisting  of  an  aggregate  of  micropolar  regions  which  satisfy  Assumptions  (i) 
and  (ii).  In  this  depiction,  there  are  Nr  =  4  regions  and  a  total  of  N  =  NrNq  Mg  cations  in  the  material. 


3.1.1.  Thermodynamics  of  Micropolar  Regions 

We  consider  first  the  behavior  of  the  micropolar  regions  which  is  determined  by  the  density  of  Mg  cations  which 
fill  B'  sites.  Under  the  assumption  that  the  regions  consists  of  a  lattice  of  cells  as  depicted  in  Figure  2,  the  dipole 
moment  p  for  the  lattice  is  given  by 

N0 

P  =  PoJ2si 

i—1 

where  po  is  the  dipole  moment  for  a  single  cell  and  No  denotes  the  number  of  Mg  cations  in  the  region.  For  PMN,  only 
cells  containing  a  Mg  cation  possess  an  associated  spin  s*  =  ±1  where  the  sign  depends  upon  the  cell’s  orientation. 
A  correlation  between  cells  is  provided  by  the  average 


.  N0 

R=w0i- 


or  equivalently,  R  = 

As  noted  in  Assumptions  (i)  and  (ii),  we  assume  that  the  regions  all  have  the  same  number  No  of  Mg  cations  but 
differing  densities  and  hence  volumes.  Within  a  region,  we  let  N+  and  N _  respectively  denote  the  number  of  positive 
and  negative  cells.  By  noting  that  N0  =  N+  +  N-,  the  cell  numbers  can  be  expressed  in  terms  of  the  correlation 
through  the  relations 

N+  =  ^(l+R)  ,  JV_  =  ^(1-R).  (6) 

The  density  X  is  defined  as  the  number  of  neighboring  B'  sites  that  contain  a  Mg  cation.  As  indicated  from 
Figure  2,  X  can  range  from  1  to  12.  Furthermore  X  is  assumed  to  be  uniform  throughout  the  region  in  accordance 
with  Assumption  (ii). 

As  will  be  detailed  later,  the  microregions  have  a  distribution  of  Curie  temperatures  due  to  the  differing  cation 
densities  in  the  regions.  To  determine  these  densities,  we  note  that  B'  sites  have  a  ratio  of  2/3  Mg  cations  to  1/3 
Nb  ions.  For  a  Mg  populated  B'  site,  each  of  the  12  adjacent  B'  sites  has  a  65.4%  probability  of  being  populated  by 
a  Mg  cation  and  34.6%  probability  of  being  populated  by  Nb.  The  probability  that  the  site  is  surrounded  by  X  Mg 
cations  is  computed  by  assuming  a  binomial  distribution 


B(X)  = 


12! 


X\{\2-X)\ 


(0.654)x  (0.346) 


12  X 


To  obtain  a  continuous  distribution  appropriate  for  subsequent  integration,  this  was  then  fit  with  the  binomial 
distribution 


npO  = 


i 


bVto 


,-(X-af/2b2 


(7) 
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having  a  mean  a  =  7.85  and  standard  deviation  b  =  1.65.  Finally,  the  distribution  is  normalized  to  the  interval  [0, 1] 
through  the  introduction  of  the  variables 


(8) 


to  obtain 

n«)  =  (9) 


The  distributions  (7)  and  (9)  will  be  employed  in  subsequent  discussion  to  compute  ferroelectric  properties  of  the 
material  which  depend  upon  cation  densities  in  the  micropolar  regions. 

To  compute  the  internal  energy  due  to  spin  conversion,  we  consider  first  the  energy  required  to  convert  an 
individual  spin  given  a  density  X  of  neighboring  Mg  cations.  Letting  $o  denote  the  energy  required  to  convert  the 
spin  at  a  single  site  when  the  system  is  completely  ordered  ( R  =  ±1),  the  energy  required  to  convert  a  positive  spin 
to  negative  (and  conversely)  is  given  by 


_  XN+ 

— - 

12  No  0 


<I> 


X  N-  ^ 

— - $0  ■ 

12  No  ° 


(10) 


We  point  out  that  these  energy  expression  are  derived  under  the  assumption  that  dipoles  interact  only  with  adjacent 
neighbors. 

The  change  in  the  internal  energy  U  due  to  spin  conversion  can  then  be  expressed  as 


dU  =  +  $_^+diV+  .  (11) 

By  utilizing  the  energy  relations  (10)  and  noting  from  (6)  that  dN+  =  ^ dR  =  — dN _,  the  expression  (11)  can  be 
integrated  to  obtain  the  relation 

U  =  ^^M(i  _  B2)  +  u0  , 

where  Uo  denotes  the  energy  for  the  completely  ordered  state,  for  the  internal  energy. 

For  stress  free  conditions,  the  total  Gibbs  free  energy  T  for  the  microregion  is  related  to  the  internal  energy  by 


r  =  u  -  STV0 


where  S  denotes  the  entropy  and  Vo  is  the  volume  of  the  microregion.  From  statistical  mechanics,  the  entropy  per 
unit  volume  for  the  system  can 


_  ksNo 

- 


ln(2) 


l±^ln(l  +  R) 


yw-*) 


which  yields  the  free  energy 

r  =  _  tf)  +  k«N»T  Lin  ( 

4  v  ;  2  [  \l-Rj 


ln(l  -  R 2) 


+  Uq  +  ksX o  ln(2) 


for  the  region.  We  note  that  the  expression  for  T  incorporates  the  ordering  of  neighboring  moments  reflected  in  the 
first  term  and  the  tendency  for  disorder  due  to  thermal  agitation  as  reflected  in  the  second  term. 

The  electric  field  Eioc  in  the  polar  region  is  determined  by 

PoNqEioc  =  ^  =  -  +  Tfcg arctanh(i?) . 

In  the  absence  of  a  local  field,  the  equation  can  be  expressed  as 

T 

R  =  — — arctanh(i?)  (12) 

-*■  cur 


where 


=  C^o 

2k  b 


(13) 


7 


denotes  the  local  Curie  temperature  for  that  region.  To  motivate  the  definition  (13)  and  this  terminology,  it  is  noted 
that  (12)  has  a  single  solution  for  T  >  Tcur  and  three  solutions  comprised  of  two  stable  and  one  unstable  solution 
for  temperatures  T  <  Tcur.  Hence  a  pitchfork  bifurcation  occurs  at  T  =  Tcur  which  is  consistent  with  a  second-order 
phase  transition. 

The  expression  (13)  quantifies  the  Curie  temperature  for  a  single  micropolar  region.  The  macroscopic  Curie 
temperature  for  the  material  is  specified  as  the  average  Curie  temperature  for  the  regions  or  the  Curie  temperature 
of  the  mean, 


T 


cur  — 


C$0 
2k  b 


(14) 


The  behavior  of  each  polar  regions  can  then  be  expressed  as 


PqNqEioc 


$0 

~2 


arctanh(i?) 


(15) 


in  the  absence  of  interactions  between  regions. 


3.1.2.  Macroscopic  Relations 


The  previous  discussion  focuses  on  the  thermodynamics  within  a  micropolar  region  but  does  not  address  interac¬ 
tions  between  regions.  As  detailed  in  [18],  however,  such  interactions  can  significantly  affect  the  effective  field  present 
within  the  material,  and  we  focus  here  on  the  development  of  commensurate  constitutive  relations  for  the  aggregate 
material  which  incorporates  these  effects. 

Following  the  development  in  [13],  we  propose  the  constitutive  relation 


|E(T)|  = 


$d3NR 
2  PSV 


~iVfi|P| 

VPs 


+ 


NR(Tf)  T 
V  Tf 


arctanh 


(16) 


To  simplify  subsequent  discussion,  we  will  consider  the  polarization  and  field  to  be  scalar- valued.  The  extensions  to 
the  vector  case  are  direct. 

We  now  employ  the  theory  developed  in  Section  3.1.1  for  the  thermodynamics  of  the  micropolar  regions  to  express 
Ps  and  Nr  in  terms  of  the  operating  temperature  and  coefficients  which  must  be  estimated  for  a  specific  material. 

The  saturation  polarization  occurs  when  the  dipole  moments  in  all  the  microregions  are  aligned  with  the  field.  The 
degree  of  alignment  for  microregions  is  quantified  by  the  correlation  R(T,X )  specified  by  (15)  while  the  magnitude 
is  dependent  upon  the  number  n(X)  of  Mg  cations  having  X  nearest  neighbors.  Recalling  that  po  denotes  the  dipole 
moment  of  a  single  cell,  the  saturation  polarization  can  be  expressed  as 


Pt 


(- T)  =  Ps  f 
Jo 


\R(T,  <)|n«)dc. 


(17) 


The  distribution  n  is  computed  using  (9)  and  the  parameter  Ps  =  NPo/V  is  estimated  either  directly  from  the  data 
or  through  a  least  squares  fit  to  the  data.  The  correlation  R  is  computed  through  iterative  solution  of  the  relation 


m0o 


arctanh(i?) 


(18) 


where  /3  =  2p£N°  denotes  the  average  degree  to  which  local  fields  turn  the  polarization  and  hence  change  the 
correlation;  see  (15).  We  note  that  at  low  drive  levels,  this  contribution  is  small  and  one  can  employ  j3  =  0.  For  high 
drive  levels,  however,  Ps  is  more  accurately  approximated  by  treating  /3  as  a  parameter  to  be  estimated  and  solving 
(18)  with  Eioc  —  Ernax. 

The  number  of  micropolar  regions  Nr,  at  a  fixed  temperature  T,  is  computed  by  determining  those  regions  with 
nonzero  correlation  R,  or  from  (15),  those  regions  in  which  T  <  Tcur.  Since  the  number  of  regions  with  density  X  is 
given  by  n(X)/N0,  the  total  number  of  regions  in  the  volume  is 


Nr(T) 

V 


n  r1 

V  No  J^T/Tcur 


i2n(c)dc- 


(19) 
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The  expressions  (17)  and  (19)  for  Ps  and  Nr  are  employed  in  (16)  to  determine  the  anhysteretic  component  of 
the  constitutive  model  (e.g.,  compare  to  (1)  and  (2)).  The  hysteresis  due  to  domain  wall  effects  is  addressed  in  the 
next  section.  The  combination  of  these  components  in  the  manner  summarized  in  Section  2  provides  the  complete 
constitutive  model  for  the  material. 


3.2.  Domain  Wall  Losses 


The  anhysteretic  model  (16)  quantifies  the  polarization  which  would  occur  in  the  absence  of  domain  wall  pinning. 
For  certain  parameter  choices,  this  model  produces  a  form  of  hysteresis  similar  to  that  observed  in  a  single  crystal 
which  undergoes  domain  switching.  However,  the  transition  from  remanence  to  the  coercive  point  is  steeper  than  that 
observed  in  most  operating  regimes.  Furthermore,  the  parameters  obtained  at  temperatures  rising  to  the  material’s 
freezing  point  typically  yields  a  single- valued  anhysteretic  curve  whereas  hysteresis  is  still  present  in  the  E-P  relation. 
These  phenomena  are  attributed  to  energy  losses  due  to  domain  wall  pinning. 

As  summarized  in  Section  2,  the  quantification  of  the  irreversible  and  reversible  polarizations,  respectively,  requires 
the  parameters  k  and  c.  We  focus  here  on  the  quantification  of  the  temperature-dependence  in  k  since  it  directly 
affects  the  decrease  in  coercivity  observed  as  temperatures  increase  to  the  freezing  point.  The  parameter  c  is  taken 
as  constant  since  its  variation  with  temperature  appears  to  be  secondary.  If  further  refinement  is  required,  analysis 
of  the  type  developed  for  k  can  be  extended  to  include  the  reversible  effects  quantified  by  c. 

For  fixed  temperature,  it  is  noted  in  Section  2  that  through  the  definition  k  =  —  ,  the  parameter  k  depends 

on  the  180°  domain  wall  energy  (£„)  required  to  move  between  potential  wells.  The  temperature  dependence  in  k  is 
attributed  to  the  temperature-dependence  of  the  potential  wells  and  hence  {£„)■  To  quantify  the  reduction  in  energy 
required  to  move  between  wells  at  increasing  temperatures,  we  employ  the  asymptotic  relation  k  «  Ec  and  focus  on 
Ec. 

The  influence  of  temperature  on  the  coercive  field  behavior  is  predicted  using  energy  arguments  drawn  from  the 
theory  of  dislocation  plasticity  [14]  and  is  similar  to  the  approach  employed  by  Chen  and  Wang  [4].  In  this  theory, 
energy  barriers  AG  due  to  local  resistive  fields  E  are  considered  as  the  mechanism  which  inhibits  domain  wall 
movement  and  hence  prohibits  domain  wall  movement  and  hence  provides  potential  wells.  Using  energy  arguments 
analogous  to  those  for  slip  dislocations,  the  application  of  Boltzmann  theory  yields  the  modified  Arrhenius  equation 


dP 

dt 


=  AP0 


Aa 

X 


Mo  exp 


-AG 

kRT 


(20) 


for  the  quasistatic  polarization  rate.  Here  Mo,  Aa  and  b  respectively  denote  the  frequency  at  which  the  domain  wall 
vibrates,  the  distance  traveled  by  the  wall,  and  the  distance  between  walls. 

To  specify  the  energy  barrier  AG,  we  consider  the  approximation 

AG  =  V*P*(E  -  E) 

employed  by  Chen  and  Wang  [4]  where  V*  denotes  the  volume  of  a  short  range  barrier  and  P*  is  the  associated 
change  in  polarization.  For  a  fixed  frequency,  this  yields 


dPdE 
dE  dt 


■  pi  exp 


-V*P*{E  -  E) 


where  pi  =  APoAap,o/b.  With  the  definition  p2  =  it  follows  that 

E  =  T-^—  in  +  £• 

Y*P*  X  dE) 

Consideration  of  the  coercive  field  E  =  Ec  then  yields 

ks 


Er  =  T- 


V*P ’ 


■\n[p2Xc{T)]+  E 


(21) 


where  Xc{T)  denotes  the  differential  susceptibility  at  Ec.  The  susceptibility  will  be  temperature  dependent  since 

Ps(T) 


Xc{T)  »  Xan{T) 


3  a{T) 
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as  indicated  in  the  asymptotic  relations  summarized  in  Section  2.2.  The  contribution  In \p2Xc{T)\  is  dominated  by 
T,  however,  so  to  a  first  approximation,  we  can  represent  the  temperature  dependence  of  the  coercive  field  by 

Ec  =  k{Tf  -  T)  (22) 

where  Tf  is  the  freezing  temperature  where  remanence  disappears.  The  linear  relation  (22)  is  consistent  with  the 
data  plotted  in  Figure  2  of  [4]  for  a  variety  of  materials.  Note  that  the  parameter  k  is  material  dependent. 

Finally,  from  the  asymptotic  relation  k  &  Ec  for  soft  materials,  we  consider  the  relation 

k{T)  =  k(Tf  -  T)  (23) 

for  the  pinning  constant  at  temperatures  T  <Tf  where  hysteresis  is  present.  For  temperatures  T  >  Tf,  the  material 
exhibits  no  hysteresis  so  k(T)  =  0  and  the  irreversible  relation 


Pirr(T)  =  Pan(T)  -  5e0k(T)-^(T) 


(24) 


reduces  to  the  anhysteretic  model. 


4.  Model  Validation 

The  modeled  derived  in  Section  3  quantifies  the  polarization  as  a  function  of  both  temperature  T  and  input  field 
E.  To  illustrate  the  predictive  capabilities  of  the  model,  we  consider  the  characterization  of  PMN-PT-BT  for  input 
fields  having  a  peak  value  of  1  MV/m  and  temperatures  ranging  from  T  =  263°  K  to  T  =  313°  K.  At  the  lower 
temperature,  the  material  exhibits  ferroelectric  behavior  and  significant  hysteresis.  As  the  temperature  increases 
through  the  freezing  temperature  (approximately  288°  K),  the  hysteresis  vanishes  and  the  material  exhibits  a  single¬ 
valued  superparaelectric  relation  between  the  applied  field  and  measured  polarization.  The  data  was  collected  from 
a  stress-free  a  —  0  sample  and  the  input  field  was  cycled  at  1  Hz  to  maintain  quasi-static  conditions. 

The  parameter  values  a  =  1.5  x  10 7,o  =  5.1,  k  =  4.0  x  103,PS  =  0.53 ,/3  =  7  x  10~9  and  c  =  .3  were  obtained 
through  a  least  squares  fit  to  the  data.  The  freezing  temperature  was  specified  to  be  Tf  =  288°  K  and  the  Curie 
temperature  was  taken  to  be  Tcur  =  313°  K.  We  reiterate  that  the  temperature-dependence  in  the  parameters 
a(T) ,  a(T) ,  k(T)  and  PS{T)  is  incorporated  solely  through  the  relations  (17)  and  the  same  parameter  set  was  employed 
through  the  full  range  of  temperatures  and  and  field  inputs. 

The  model  fit  obtained  with  these  parameters  is  compared  with  the  measured  data  in  Figure  5.  It  is  observed 
that  the  model  very  accurately  quantifies  the  initial  hysteresis  at  T  =  263°  K,  the  decreased  hysteresis  at  T  =  283° 
K,  and  the  anhysteretic  behavior  at  higher  temperatures.  It  also  quantifies  the  decrease  in  the  maximum  measured 
polarization  which  occurs  as  temperatures  increase.  We  point  out  the  model  parameters  were  chosen  to  optimize 
the  model  predictions  throughout  the  full  temperature  range  and  better  can  be  obtained  for  specific  temperature 
intervals  if  desired.  For  example,  the  parameter  choices  a  =  1.4  x  107,a  =  4.1,  A;  =  5.0  x  103,PS  =  0.53  and  c  =  .5 
yield  a  highly  accurate  model  fit  to  the  hysteresis  at  T  =  263°  K  but  with  less  accuracy  at  the  higher  temperatures. 

5.  Concluding  Remarks 

This  paper  addresses  the  modeling  of  temperature-dependent  constitutive  nonlinearities  and  hysteresis  in  relaxor 
ferroelectrics  materials.  The  model  is  based  on  the  assumption  that  the  materials  are  comprised  of  micropolar  regions 
having  a  broad  spectrum  of  Curie  temperature  due  to  their  chemical  heterogeneity.  Consideration  of  the  free  energy 
in  the  presence  of  variations  in  cation  distributions  yields  an  Ising  spin  model  representing  the  thermodynamics 
of  the  microregions.  Based  on  a  “random  layer”  model  for  B-site  orientations,  this  then  provided  a  mechanism 
for  computing  the  saturation  polarization  and  distribution  of  micropolar  regions  as  a  function  of  temperature.  A 
nonlinear  anhysteretic  model  for  the  aggregate  material  was  obtained  by  combining  the  micropolar  model  with 
macroscopic  constitutive  relations.  Finally,  hysteresis  was  incorporated  through  the  quantification  of  energy  required 
to  bend  and  translate  domain  walls  pinned  at  inclusions  inherent  to  the  material. 

The  predictive  capabilities  of  the  model  were  illustrated  through  comparison  with  PMN-PT-BT  data  generated 
through  the  input  of  a  1  Hz  field  having  a  peak  magnitude  of  1  MV/m  at  temperatures  ranging  from  263°  K  to 
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Figure  5.  PMN-PT-BT  data  ( - )  and  the  domain  wall  model  ( - )  for  the  temperature  range  T  =  263°K  to 

T  =  313°K.  In  each  case,  the  field  E  (MV/m)  is  plotted  versus  polarization  P  (C/m2). 


313°  K.  At  the  lower  temperatures,  the  model  quantifies  the  inherent  hysteresis.  It  also  quantifies  the  transition  to 
anhysteretic  behavior  and  the  reduction  in  saturation  polarization  which  occur  as  temperatures  increase.  The  full 
temperature  and  field  dependence  is  incorporated  in  the  model  and  the  differing  model  responses  were  obtained  with 
one  set  of  parameters.  The  sole  input  to  the  model  is  the  input  field  and  the  operating  temperature.  The  small 
number  of  parameters  (6)  and  their  asymptotic  relation  to  measurable  quantifies  (e.g.,  coercive  field  Ec.  slope  at  field 
reversal)  facilitates  the  construction  of  the  model  as  well  as  the  updating  of  parameters  to  accommodate  changing 
operating  conditions. 
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